1 Introduction

Functional principal component analysis is a common tool to explore the directions of variation in functional data. Electroencephalography (EEG) studies create data with a complex structure including regional, functional, and longitudinal dimensions and without proper tools to study this type of data, standard analysis involves averaging over one or more of these dimensions. While there are many recent developments in FPCA, many of them consider either longitudinally or spatially repeated functional data. The Hybrid Principal Component Analysis (HPCA) decomposition proposed by Scheffler et al. (2020) uses weak separability to combine vector and functional principal components analysis, building off of Product FPCA of Chen et al. (2016) which had a product one-dimensional eigenfunctions.

2 Hybrid Principal Component Analysis

Scheffler et al. (2020) proposed Hybrid Principal Component Analysis (HCPA), a method for analyzing multi-dimensional EEG data utilizing weak separability.

Let \(Y_{di}(r, \omega, s)\) denote the log principal power, which is the functional observation that we consider, for subject \(i, i = 1, \ldots, n_d\), from group \(d, d = 1, \ldots, D\), in region \(r, r = 1, \ldots, R\), at frequency \(\omega, \omega \in \Omega\), and segment \(s, s \in \mathcal{S}\).

\[\begin{align*} Y_{di}(r, \omega, s) = \mu(\omega, s) + \eta_d(r, \omega, s) + Z_{di}(r, \omega, s) + \epsilon_{di}(r, \omega, s) \end{align*}\]

Using weak separability of the covariances, the common eigenfunctions and eigenvectors are found using the marginal covariances. The functional and longitudinal marginal covariance surfaces are defined as \[\begin{align*} \Sigma_{d, \Omega}(\omega, \omega') &= \sum_r \int_\mathcal{S} \text{cov}\{Z_{di}(r, \omega', s), Z_{di}(r, \omega, s)\}ds = \sum_{\ell=1}^\infty \tau_{d\ell, \Omega} \phi_{d\ell}(\omega) \phi_{d\ell}(\omega'),\\ \Sigma_{d, \mathcal{S}}(s, s') &= \sum_r \int_\Omega \text{cov}\{Z_{di}(r, \omega, s), Z_{di}(r, \omega, s')\}d\omega = \sum_{m=1}^\infty \tau_{dm, \mathcal{S}}\psi_{dm}(s)\psi_{dm}(s'), \end{align*}\] and the regional marginal covariance matrix with \((r, r')\)-th element is \[\begin{align*} (\Sigma_{d, \mathcal{R}})_{r, r'} &= \int_\mathcal{S}\int_\Omega \text{cov}\{Z_{di}(r, \omega, s), Z_{di}(r', \omega, s)\}d\omega ds = \sum_{k=1}^R \tau_{dk, \mathcal{R}}\text{v}_{dk}(r)\text{v}_{dk}(r') \end{align*}\] where \(\phi_{d\ell}(\omega)\) and \(\psi_{dm}(s)\) are the common eigenfunctions of the functional and longitudinal marginal covariance surfaces and \(v_{dk}(r)\) are the common eigenvectors for the regional marginal covariance matrix, each with their respective eigenvalues \(\tau_{d\ell, \Omega}\), \(\tau_{dm}(s)\) and \(\tau_{dk, \mathcal{R}}\).

Now the HPCA decomposition of \(Y_{di}(r, \omega, s)\) is given as \[\begin{align*} Y_{di}(r, \omega, s) &= \mu(\omega, s) + \eta_d(r, \omega, s) + \sum_{k=1}^R\sum_{\ell=1}^\infty\sum_{m=1}^\infty \xi_{di, k\ell m} \text{v}_{dk}(r) \phi_{d\ell}(\omega)\psi_{dm}(s) + \epsilon_{di}(r, \omega, s), \end{align*}\] where the subject-specific scores \(\xi_{di, k\ell m}\) are defined through the projection \(\langle Z_{di}(r, \omega, s), v_{dk}(r)\phi_{d\ell}(\omega)\psi_{dm}(s) \rangle = \sum_{r=1}^R\int\int Z_{di}(r,\omega, s)v_{dk}(r)\phi_{d\ell}(\omega)\psi_{dm}(s) d\omega ds.\) The decomposition of the total covariance is \[\begin{align*} \Sigma_{d, T}\{(r, \omega, s), (r', \omega', s')\}&= \text{cov}\{Z_{di}(r, \omega, s), Z_{di}(r', \omega', s')\} + \sigma^2_d\delta\{(r, \omega, s), (r', \omega', s')\}\\ &= \sum_{k=1}^R\sum_{\ell=1}^\infty\sum_{m=1}^\infty \tau_{d, k\ell m} \text{v}_{dk}(r)\phi_{d\ell}(\omega)\psi_{dm}(s) \text{v}_{dk}(r')\phi_{d\ell}(\omega')\psi_{dm}(s') + \sigma^2_d\delta\{(r, \omega, s), (r', \omega', s')\} \end{align*}\] where \(\tau_{d, k\ell m} = \text{var}(\xi_{di, k\ell m})\) and \(\delta\{(r, \omega, s), (r', \omega', s')\}\) denotes the indicator for \(\{(r, \omega, s) = (r', \omega', s')\}\).

The three-dimensinoal HPCA can be reduced to a two-dimensional HPCA when the longitudinal dimension is not of interest. The decomposition of \(Y_{di}(r, \omega)\) then becomes \[\begin{align*} Y_{di}(r, \omega) &= \mu(\omega) + \eta_d(r, \omega) + Z_{di}(r, \omega) + \epsilon_{di}(r, \omega) \\ &= \mu(\omega) + \eta_d(r, \omega)+\sum_{k=1}^R\sum_{\ell=1}^\infty \xi_{di, k\ell} \text{v}_{dk}(r) \phi_{d\ell}(\omega)+ \epsilon_{di}(r, \omega). \end{align*}\] This will be of important in our data application since the data I have avaiable has already been averaged over trials.

Without repeating all of the details in the article, a summary of the HPCA estimation procedure’s steps are given below.

Algorithm: HPCA estimation procedure
1. Estimation of fixed effects:
2. Estimation of marginal covariances and measurement error variance
3. Estimation of marginal eigencomponents
4. Estimation of variance components and subject-specific scores via linear mixed effects models

The triple index \(k\ell m\) in HPCA truncated at $K, L, $ and \(M\) is replaced with a single index \(g = (k - 1) + K(\ell - 1) + KL(m - 1)\) and the model becomes \[\begin{align*} Y_{di}(r, \omega, s) = \mu(\omega, s) + \eta_d(r, \omega, s) + \sum_{g = 1}^G \zeta_{dig}\varphi_{dg}(r, \omega, s) + \epsilon_{di}(r, \omega, s), \end{align*}\] where \(\varphi(r, \omega, s) = v_{dk}(r)\phi_{d\ell}(\omega)\psi_{dm}(s)\), \(\zeta_{dig} = \langle Z_{di}(r, \omega, s), \varphi(r, \omega, s) \rangle\), \(\kappa_{dg} = \text{var}(\zeta_{dig})\), and \(G = KLM\). Now, with a vectorized version of \(Y_{di(r, \omega, s)}\) denoted by \(\boldsymbol{Y}_{di}\) and similar vectorized versions of \(\varphi_{dg}(r, \omega, s)\), \(\mu(t)\), \(\eta_d(r, t)\), and \(\epsilon_{di}(r, \omega, s)\), the linear mixed effects model is \[\begin{align*} \boldsymbol{Y}_{di} &= \boldsymbol{\mu}_i + \boldsymbol{\eta}_{di} + \boldsymbol{Z}_{di} + \boldsymbol{\epsilon}_{di} =\boldsymbol{\mu}_i + \boldsymbol{\eta}_{di} + \sum_{g=1}^{G'} \zeta_{dig} \boldsymbol{\varphi}_{dig} + \boldsymbol{\epsilon}_{di} \end{align*}\] for \(i = 1, \ldots, n_d\).

3 Group-level inference

When constructing these EEG studies, the interest lies in the contrast between conditions between groups so the group-level inference we want to conduct is testing the null hypothesis \(H_0: \eta_d(r, \omega, s) = \eta(r, \omega, s)\) for \(d = 1, \ldots, D\). Scheffler et al. (2020) proposed a parametric bootstrap procedure that generates outcomes based on the model components estimated using the HPCA estimation procedure. Generate outcomes under the null in region \(r\) as \[\begin{align*} Y_{di}^b(r, \omega, s) &= \widehat{\mu}(\omega, s) + \widehat{\eta}(r, \omega, s) + \sum_{g=1}^{G'} \zeta_{dig}^b\widehat{\varphi}_{dig}(r, \omega, s) + \epsilon_{di}^b(r, \omega, s) \end{align*}\] and in the other regions as \[\begin{align*} Y_{di}^b(r, \omega, s) &= \widehat{\mu}(\omega, s) + \widehat{\eta}_d(r, \omega, s) + \sum_{g=1}^{G'} \zeta_{dig}^b\widehat{\varphi}_{dig}(r, \omega, s) + \epsilon_{di}^b(r, \omega, s) \end{align*}\] where \(\zeta_{di}^b \sim N(0, \widehat{\kappa}_{dg})\), \(\epsilon_{di}^b(r, \omega, s) \sim N(0, \widehat{\sigma}^2_d)\), and \(\widehat{\eta}(r, \omega, s)\) is the point-wise average of the group-region shift estimates. Then the proposed test statistic is \(T_r = [\sum_{d=1}^D \int\int \{\widehat{\eta}_d(r, \omega, s) - \widehat{\eta}(r, \omega, s)\}]^{1/2}\). The statistic is calculated for each bootstrap simulation and the p-value is calculated as $_{b = 1}^B (T_R^b > T_R).

4 Data Application

The first application of pERP-RED examines EEG data from a study of the neural mechanisms underlying language impairment in children with Autism Spectrum Disorder (ASD) detailed in DiStefano et al. (2019). This group has been difficult to study using traditional paradigms that require following instructions or providing behavioral responses. Thirty-one children, aged 5-11 years old, who participated in the below described paradigms and provided sufficient high quality EEG data were considered: 14 typically developing (TD), 10 verbal ASD (vASD), and 7 minimally verbal ASD (mvASD). The three groups were age-matched. In the ASD participants, diagnoses had been made prior to enrollment, through clinical diagnosis by independent clinical psychologists, child psychiatrists, and/or developmental pediatricians. These diagnoses were confirmed by the research team using the Autism Diagnostic Observation Schedule (ADOS) and Social Communication Questionnaire (SCQ).

Two picture-word matching paradigms were used, one audio and one visual. The word stimuli included 60 basic nouns taken from the MacArthur-Bates Communicative Development Inventories-2nd edition. Examples of words included animals such as a bird or a dog and inanimate objects such as a doll or a bike. In the audio paradigm, a picture of the word would appear on a white background and the child would hear a spoken word that was either the same (match condition) or a word neither semantically nor phonologically matched the word image (mismatch condition). In both conditions, the picture image appeared for 2000ms, where the auditory stimulus was played after 1000ms after the picture image was shown (see Figure 4.1). Each word image pair appeared twice—once in the matched and once in the mismatched condition. No behavioral response was required. Trials were presented in four blocks of 30 trials each, totaling 6 minutes.

The visual paradigm used the same nouns and same number of trials, but in each an image of the word appeared on screen after 1000ms rather than an audio recording (see Figure 4.1). As in the auditory paradigm, the two conditions are ‘’match’’ or ‘’mismatch’’. In both paradigms, the trials were video recorded in order to remove trials in which the participants were not looking at the screen.

In order to use the HPCA decomposition, we will look at the contrast of the match and mismatch conditions of the audio paradigm. This paradigm was chosen rather than the visual paradigm because it has been studied more extensively in the DiStefano et al. (2019) and Campos et al. (2020) publications. As this is an event-related study, the functional domain is time rather than frequency as in Scheffler et al. (2020). The main findings from previous works include a deeper negativity for mismatch than match trials, which they labeled as an N4 and late negative component (LNC). It was also found that there was more heterogeneity in the ASD groups than the TD group.

In both paradigms, the image appears for the full 2000ms. (a) In the audio paradigm, a spoken word is heard 1000ms after the onset of the image. (b) In the visual paradigm, the word image is seen 1000ms after the onset of the image.In both paradigms, the image appears for the full 2000ms. (a) In the audio paradigm, a spoken word is heard 1000ms after the onset of the image. (b) In the visual paradigm, the word image is seen 1000ms after the onset of the image.

Figure 4.1: In both paradigms, the image appears for the full 2000ms. (a) In the audio paradigm, a spoken word is heard 1000ms after the onset of the image. (b) In the visual paradigm, the word image is seen 1000ms after the onset of the image.

The estimated mean function \(\widehat{\mu}(t)\) plotted in Figure 4.2 shows that in general the electrodes do not have a common trend, since the overal mean function is nearly identical to 0. Figure 4.3 depicts the estimated \(\eta_d(r, t)\) functions for each region and group. From this, it is possible to see where there may be differences when performing group-level inference. The plots for the C3, P10, and P3 electrodes show some group-level deviation. Table 4.1 shows the proportion of variation explained by each of the eigencomponents for the eigencomponents that explain at least 90% variation in all three diagnostic groups.

The estimated mean function (blue) is identically zero.

Figure 4.2: The estimated mean function (blue) is identically zero.

The group-region deviations from the mean are the subject of inference in this study. There is evidence that the $\eta$ functions may differ for the P10, C4, and CZ electrodes across groups.

Figure 4.3: The group-region deviations from the mean are the subject of inference in this study. There is evidence that the \(\eta\) functions may differ for the P10, C4, and CZ electrodes across groups.

Table 4.1: FVE of the marginal covariances
TD
vASD
mvASD
T R T R T R
0.3081450 0.4294286 0.2936115 0.4810002 0.4373479 0.5234623
0.2005062 0.2546679 0.2128275 0.1841818 0.1911425 0.2075404
0.1629697 0.1569455 0.2055434 0.1437301 0.1121524 0.1513102
0.1234887 0.0877757 0.0864545 0.0620644 0.0819160 0.0554136
0.0831330 0.0803797 0.0555993 0.0627756
0.0372993 0.0431772 0.0335934

The first two estimated eigenfunctions are plotted in Figure 4.4. The first leading eigenfunction shows the most variation from 200 to 750 ms in the TD and vASD groups and from 500 to 900 ms in the mvASD group. These are the two ranges of time that were of interest in the DiStefano et al. (2019) paper. This is evidence of the heterogeneity in the semantic integration response. The first two estimated eigenvectors are plotted in Figure 4.5. The first leading eigenvector shows a contrast between the front and the back of the head for all diagnostic groups while the second leading eigenvector shows a contrast between the left and right sides of the head in all groups.

Estimated first and second leading functional eigenfunctions, $\phi_{d1}(t)$, $\phi_{d2}(t)$.Estimated first and second leading functional eigenfunctions, $\phi_{d1}(t)$, $\phi_{d2}(t)$.

Figure 4.4: Estimated first and second leading functional eigenfunctions, \(\phi_{d1}(t)\), \(\phi_{d2}(t)\). 

electrodeLocs <- read_csv(here("data/electrodeLocs81.csv"), 
                          col_types = cols())

theme_topo <- function(base_size = 12) {
  theme_bw(base_size = base_size) %+replace%
    theme(
      rect = element_blank(),
      line = element_blank(),
      axis.text = element_blank(),
      axis.title = element_blank(),
      legend.position = "bottom"
    )
}

circleFun <- function(center = c(0, 0), diameter = 1, npoints = 100) {
  r = diameter / 2
  tt <- seq(0, 2 * pi, length.out = npoints)
  xx <- center[1] + r * cos(tt)
  yy <- center[2] + r * sin(tt)
  return(data.frame(x = xx, y = yy))
}

headShape <- circleFun(c(0, 0), round(max(electrodeLocs$x)), npoints = 100)
nose <- data.frame(x = c(-0.075, 0, .075), y = c(.495, .575, .495))

jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", 
                                 "#7FFF7F", "yellow", "#FF7F00", "red", 
                                 "#7F0000"))

limits <- map_dfr(
  data.frame(
    mvASD = -soundcontrast_HPCA$reg$mvASD$eigendecomposition$vectors[, 1],
    vASD  = -soundcontrast_HPCA$reg$vASD$eigendecomposition$vectors[, 1],
    TD    =  soundcontrast_HPCA$reg$TD$eigendecomposition$vectors[, 1]
  ),
  ~ data.frame(.x) %>%
    setNames("v") %>% 
    mutate(region = regional),
  .id = "Group"
) %>%
  mutate(Group = fct_relevel(Group, "TD", "vASD", "mvASD")) %>% 
  ungroup() %>% 
  summarise("max" = max(v), 
            "min" = min(v))

scores_long <- map_dfr(
  data.frame(
    mvASD = -soundcontrast_HPCA$reg$mvASD$eigendecomposition$vectors[, 1],
    vASD  = -soundcontrast_HPCA$reg$vASD$eigendecomposition$vectors[, 1],
    TD    =  soundcontrast_HPCA$reg$TD$eigendecomposition$vectors[, 1]
  ),
  ~ data.frame(.x) %>%
    setNames("v") %>% 
    mutate(electrode = regional),
  .id = "Group"
) %>%
  mutate(Group = fct_relevel(Group, "TD", "vASD", "mvASD")) %>% 
  ungroup() %>%
  left_join(electrodeLocs, by = "electrode")

ggplot(headShape, aes(x, y)) +
  geom_path(size = 1.5) +
  geom_point(data = scores_long,
             aes(x, y, colour = v),
             size = 7) +
  geom_text(data = scores_long, aes(x, y, label = electrode), size = 3) + 
  scale_colour_hp(
    option = "LunaLovegood",
    guide = "colourbar",
    oob = squish, 
    limits = c(limits$min, limits$max)
  ) +
  geom_line(data = nose, aes(x, y, z = NULL), size = 1.5) +
  theme_topo() +
  facet_wrap(~ Group, ncol = 5) +
  coord_equal()  + 
  labs(
    title = "First leading eigenvector",
    color = TeX(glue("$v_{1}(r)$"))
  )


limits <- map_dfr(
  data.frame(
    mvASD =  soundcontrast_HPCA$reg$mvASD$eigendecomposition$vectors[, 2],
    vASD  = -soundcontrast_HPCA$reg$vASD$eigendecomposition$vectors[, 2],
    TD    =  soundcontrast_HPCA$reg$TD$eigendecomposition$vectors[, 2]
  ),
  ~ data.frame(.x) %>%
    setNames("v") %>% 
    mutate(region = regional),
  .id = "Group"
) %>%
  mutate(Group = fct_relevel(Group, "TD", "vASD", "mvASD")) %>% 
  ungroup() %>% 
  summarise("max" = max(v), 
            "min" = min(v))

scores_long <- map_dfr(
  data.frame(
    mvASD =  soundcontrast_HPCA$reg$mvASD$eigendecomposition$vectors[, 2],
    vASD  = -soundcontrast_HPCA$reg$vASD$eigendecomposition$vectors[, 2],
    TD    =  soundcontrast_HPCA$reg$TD$eigendecomposition$vectors[, 2]
  ),
  ~ data.frame(.x) %>%
    setNames("v") %>% 
    mutate(electrode = regional),
  .id = "Group"
) %>%
  mutate(Group = fct_relevel(Group, "TD", "vASD", "mvASD")) %>% 
  ungroup() %>%
  left_join(electrodeLocs, by = "electrode")

ggplot(headShape, aes(x, y)) +
  geom_path(size = 1.5) +
  geom_point(data = scores_long,
             aes(x, y, colour = v),
             size = 7) +
  geom_text(data = scores_long, aes(x, y, label = electrode), size = 3) + 
  scale_colour_hp(
    option = "LunaLovegood",
    guide = "colourbar",
    oob = squish, 
    limits = c(limits$min, limits$max)
  ) +
  geom_line(data = nose, aes(x, y, z = NULL), size = 1.5) +
  theme_topo() +
  facet_wrap(~ Group, ncol = 5) +
  coord_equal()  + 
  labs(
    title = "Second leading eigenvector",
    color = TeX(glue("$v_{2}(r)$"))
  )
Estimated first and second leading regional eigenvectors, $v_{d1}(r)$ and $v_{d}(r)$.Estimated first and second leading regional eigenvectors, $v_{d1}(r)$ and $v_{d}(r)$.

Figure 4.5: Estimated first and second leading regional eigenvectors, \(v_{d1}(r)\) and \(v_{d}(r)\).

As in the Scheffler et al. (2020) paper, we utilize the 2-dimensional HPCA bootstrap adaptation via the test statistic \(T_r = [\sum_{d=1}^D \int \{ \widehat{\eta}_d(r, t) - \widehat{\eta}(r, t)\} dt]^{1/2}\). For each scalp region, we test the null hypothesis \(H_0: \eta_d(r, t) = \eta(r, t), d = 1, 2, 3\) and we find significant differences in these regions across the full time domain: C3, F7, F9, P10, P3, P8, PZ, T10, T7, see Figure 4.6.

The differences of the estimated group-region shifts $\eta_d(r, t)$ from group-region averages.

Figure 4.6: The differences of the estimated group-region shifts \(\eta_d(r, t)\) from group-region averages.

5 Discussion

The hybrid principal component analysis decomposition is an effective tool to characterize the variation in multiple dimensions without needing to collapse across an important dimension. The estimation procedure handles large data sets efficiently and can be adapted for 2 and 3 dimensional structures. Currently, I am developing a multilevel model which incorporates the ideas of weak separability in the same ways that HPCA incorporated it. The model has the form: \[\begin{equation*} Y_{dij}(r, \omega, s) = \mu(\omega, s) + \eta_d(r, \omega, s) + Z_{di}(r, \omega, s) + W_{dij}(r, \omega, s) + \epsilon_{dij}(r, \omega, s) \end{equation*}\] and will allow researchers to characterize the variation between observations in data applications such as the one described above or a multi-visit EEG study.

6 References

Campos, E., Hazlett, C., Tan, P., Truong, H., Loo, S., DiStefano, C., Jeste, S., & Senturk, D. (2020). NeuroImage Principle ERP reduction and analysis : Estimating and using principle ERP waveforms underlying ERPs across tasks , subjects and electrodes. NeuroImage, 212(August 2019). https://doi.org/10.1016/j.neuroimage.2020.116630

Chen, K., Delicado, P., & Muller, H.-G. (2016). Modeling Function-Valued Stochastic Processes , With Applications to Fertility Dynamics. Journal of the Royal Statistical Society: Series B (Statistical Methodology, 79, 177–196.

DiStefano, C., Senturk, D., & Jeste, S. S. (2019). ERP evidence of semantic processing in children with ASD. Developmental Cognitive Neuroscience, 36, 100640. https://doi.org/10.1016/j.dcn.2019.100640

Scheffler, A., Telesca, D., Li, Q., Sugar, C. A., Distefano, C., Jeste, S., & Senturk, D. (2020). Hybrid principal components analysis for region-referenced longitudinal functional EEG data. Biostatistics, 21(1), 139–157. https://doi.org/10.1093/biostatistics/kxy034

7 Code Appendix

Code for the 2D HPCA decomposition:

# Packages and Scripts ---------------------------------------------------

library(tidyverse)  # loads many packages with many useful tools
library(here)       # package for path reproducibility
library(Matrix)     # used for sparseMatrix
library(tictoc)     # timing code
library(data.table) #
library(mgcv)       # GAMs 
library(fANCOVA)    #
library(pracma)     # pracma::trapz -- trapezoidal integration
library(rlist)      # rlist::list.rbind -- rbind a list of dataframes


# My HPCA Adaptation ------------------------------------------------------

HPCA_2d <- function(data, functional, regional, kCOV, reControl) { 
  
  tictoc::tic("Estimation")
  
  # 0. Format data and create return list ----------------------------------
  
  tictoc::tic("0. Data formatting")
  
  # format data and sort
  data <- data.table(data) 
  data <- arrange(data, Group, Subject, Electrode, Time) 
  
  # define global variables 
  id                <- unique(data$Subject)
  groups            <- unique(data$Group)
  names(groups)     <- groups
  d_tot             <- length(groups)
  f_tot             <- length(functional)
  r_tot             <- length(regional) 
  
  # create dimension indices for merging
  data <- data %>%
    mutate(
      func_index = match(Time, functional),
      reg_index  = match(Electrode, regional),
      index      = (reg_index - 1) + r_tot * (func_index - 1) + 1
    )
  
  # create list for output 
  HPCA <- matrix(list(), 6, 1) 
  names(HPCA) <- c("mu", "eta", "functional", "reg", "RE model", "data")
  
  tictoc::toc()
  
  
  # 1. Estimation of Fixed Effects -----------------------------------------
  
  tictoc::tic("1. Estimation of Fixed Effects")
  
  # a. Calculate overall mean function
  
  # calculate overall mean function within each group
  HPCA[[1]] <- 
    # calculate overall mean function within each group
    map_dfc(
      groups, 
      function(g) {
        dat <- filter(data, Group == g)
        spline_fit <- smooth.spline(dat$Time, dat$ERP)
        group_mean <- predict(spline_fit, functional)$y
      }
    ) %>% 
    # add in Time variable 
    mutate(Time = functional) %>% 
    # rearranging variables 
    select(Time, everything()) %>% 
    # calculate overall mean function as the pointwise average of the group
    # means
    mutate(overall = (TD + vASD + mvASD)/d_tot) %>% 
    pull(overall) 
  
  # center data by overall mean function 
  data <- data %>% 
    group_by(Subject, Group, Electrode) %>% 
    mutate(Overall_Mean = HPCA[[1]], 
           ERP_centered = ERP - Overall_Mean)
  
  # b. calculate group-region-specific shifts 
  
  # obtain group-region-specific shifts 
  HPCA[[2]] <- map_dfr(
    groups, 
    function(g) { 
      map_dfc(
        regional, 
        function(r) {
          dat <- filter(data, Group == g, Electrode == r)
          spline_fit <- smooth.spline(dat$Time, dat$ERP_centered)
          group_region_mean <- predict(spline_fit, functional)$y
        }
      ) %>% 
        mutate(Group = g, Time = functional) %>% 
        select(Group, Time, everything())
    }
  ) %>% 
    gather(Electrode, Eta_Mean, -c(Group, Time))
  
  # center data by group-region-specific shifts
  data <- full_join(data, 
                    HPCA[[2]], 
                    by = c("Group", "Electrode", "Time")) %>% 
    mutate(ERP_centered2 = ERP_centered - Eta_Mean) %>% 
    ungroup()
  
  tictoc::toc()
  
  
  # 2. Estimation of Marginal Covariances and Measurement Error Variance ----
  
  tictoc::tic("2. Estimation of Marginal Covariances and Measurement Error Variance")
  
  #' @title Marginal Covariance Estimator 
  #' @description function for obtaining marginal covariances across regional 
  #' or functional dimensions 
  #'
  #' @param data data.table from preceding estimation 
  #' @param g group indicator 
  #' @param marginal_domain marginal domain 
  #' @param marginal_index marginal domain index 
  #' @param fixed_index fixed domain index 
  #' @param marginal_grid number of grid points in marginal domain 
  #' @param fixed_grid number of grid poitns in fixed domain 
  #' @param smooth TRUE = covariance smoothing, FALSE otherwise 
  #' @param k knots for tensor spline basis, scalar
  #'
  #' @return
  
  marginal_covariance <- function(data, g, marginal_domain, marginal_index,
                                  fixed_index, marginal_grid, fixed_grid, 
                                  smooth, k) {
    
    # subset data by group 
    data <- filter(data, Group == g) %>% 
      ungroup()
    
    # obtain unique group IDs
    group_subjects <- unique(data$Subject) 
    
    # number of subjects in group (n_d) 
    n_d <- length(group_subjects) 
    
    # assign an index to the columns of the marginal design matrix 
    data <- mutate(
      data,
      row_index = (match(Subject, group_subjects) - 1) +
        n_d * (!!sym(fixed_index) - 1) + 1
    )
    
    # form sparse marginal design matrix 
    sparse_mat = sparseMatrix( 
      i = pull(data, row_index), 
      j = pull(data, marginal_index), 
      x = pull(data, ERP_centered2), 
      dims = c(n_d * fixed_grid, marginal_grid)
    )
    
    # form sparse marginal indicator matrix 
    sparse_mat_ind <- sparseMatrix( 
      i = pull(data, row_index), 
      j = pull(data, eval(marginal_index)), 
      x = 1, 
      dims = c(n_d * fixed_grid, marginal_grid)
    ) 
    
    # calculate covariance matrix 
    cov_mat <- as.matrix(crossprod(sparse_mat, sparse_mat)
                         / crossprod(sparse_mat_ind, sparse_mat_ind))
    
    if (smooth == TRUE) { 
      # use for marginal functional covariance
      
      # vectorize covariance for smoothing 
      cov_vec <- data.frame( 
        covariance = as.vector(cov_mat), 
        w1 = (rep(marginal_domain, marginal_grid)), 
        w2 = (rep(marginal_domain, each = marginal_grid))
      )
      
      # remove entries without any observations 
      cov_vec <- cov_vec[complete.cases(cov_vec), ]
      
      # remove diagonal entries 
      cov_vec_d <- cov_vec[which(cov_vec$w1 != cov_vec$w2), ]
      
      # smooth the pooled sample covariance 
      cov_vec_s <- gam(covariance ~ te(w1, w2, k = 10, bs = "ps"), 
                       data = cov_vec_d)
      
      # form 2D grid to predict covariance function 
      newdata <- data.frame( 
        w1 = (rep(marginal_domain, marginal_grid)), 
        w2 = (rep(marginal_domain, each = marginal_grid))
      )
      
      # estimated marginal covariance function 
      cov_mat_s <- matrix(predict(cov_vec_s, newdata = newdata), 
                          nrow = marginal_grid)
      
      # symmetrize covariance function 
      cov_mat_s <- (cov_mat_s + t(cov_mat_s)) / 2
      
      # calculate measurement error variance
      # extract diagonal entries from pooled sample covariance 
      marg_diag <- cov_vec[which(cov_vec$w1 == cov_vec$w2), ]
      
      # smooth diagonal entries 
      loess_diag <- suppressWarnings(loess.as(
        marg_diag[, 2], 
        marg_diag[, 1], 
        degree = 1, 
        criterion = "gcv", 
        user.span = NULL, 
        plot = FALSE
      ))
      
      # calculate measurement error variance  
      sigma_2 <- mean(predict(loess_diag, marginal_domain) - diag(cov_mat_s)) 
      return(list(cov_mat, cov_mat_s, sigma_2)) 
      
    } else { 
      
      # use for marginal regional covariance 
      cov_mat_s = cov_mat - diag(HPCA[[3]][[g]][[1]][[3]], r_tot)
      
      # symmetrize covariance function
      cov_mat_s = (cov_mat_s + t(cov_mat_s)) / 2
      return(list(cov_mat, cov_mat_s))
    }
  }
  
  # estimate marginal covariances for each group 
  for (g in groups) { 
    HPCA[[3]][[g]] <- matrix(list(), 2, 1)
    names(HPCA[[3]][[g]]) <- c("covariance", "eigendecomposition")
    HPCA[[4]][[g]] <- matrix(list(), 2, 1) 
    names(HPCA[[4]][[g]]) <- c("covariance", "eigendecomposition")
    
    # a. calculate functional marginal covariance 
    HPCA[[3]][[g]][[1]] <- marginal_covariance(
      data, 
      g, 
      functional, 
      "func_index",
      "reg_index", 
      f_tot, 
      r_tot, 
      TRUE, 
      kCOV
    )
    names(HPCA[[3]][[g]][[1]]) <- c("sigma_hat", "sigma_tilde", "sigma_2")
    
    # b. calculate regional marginal covariance 
    HPCA[[4]][[g]][[1]] <- marginal_covariance(
      data, 
      g, 
      regional, 
      "reg_index", 
      "func_index", 
      r_tot, 
      f_tot,
      FALSE, 
      NA
    ) 
    names(HPCA[[4]][[g]][[1]]) <- c("sigma_hat", "sigma_tilde")
  }
  
  tictoc::toc()
  
  
  # 3. Estimation of Marginal Eigencomponents ------------------------------
  
  tictoc::tic("3. Estimation of Marginal Eigencomponents")
  
  #' @title Eigenfunctions 
  #' @descriptions function for obtaining marginal eigencomponents that explain
  #' 85% FVE 
  #'
  #' @param covariance marginal covariance function
  #' @param marginal_domain marginal domain
  #' @param type TRUE for smooth covariance function, FALSE for covariance 
  #' matrix
  #'
  #' @return
  
  eigenfunctions <- function(covariance, marginal_domain, type) { 
    
    eigen_temp <- eigen(covariance, symmetric = TRUE) 
    
    # obtain positive eigenalues 
    eigen_temp$values <- eigen_temp$values[which(eigen_temp$values > 0)] 
    
    # obtain eigenvectors associated with positive eigenvalues 
    eigen_temp$vectors <- eigen_temp$vectors[, 1:length(eigen_temp$values)]
    
    if (type == TRUE) {
      for (e in 1:length(eigen_temp$values)) {
        # normalize the eigenvalues over the domain 
        eigen_temp$vectors[, e] <- eigen_temp$vectors[, e] / 
          sqrt(trapz(marginal_domain, eigen_temp$vectors[, e]^2)) 
      }
    }
    
    # select components that explain at least 85% FVE 
    K <- length(which(cumsum(eigen_temp$values) / 
                        sum(eigen_temp$values) < .9)) + 1
    
    return(list(
      "values" = eigen_temp$values, 
      "vectors" = eigen_temp$vectors, 
      "FVE" = K 
    ))
  }
  
  # estimate eigencomponents for each group 
  for (g in groups) { 
    # employ FPCA on functional marginal covariance 
    HPCA[[3]][[g]][[2]] <- eigenfunctions(HPCA[[3]][[g]][[1]][[2]], 
                                          functional, 
                                          TRUE)
    
    # employ PCA on regional marginal covariance
    HPCA[[4]][[g]][[2]] <- eigenfunctions(HPCA[[4]][[g]][[1]][[2]], 
                                          regional, 
                                          FALSE) 
  }
  
  tictoc::toc()
  
  
  # 4. Estimation of Variance Components and Subject-specific Scores -------
  
  tictoc::tic("4. Estimation of Variance Components and Subject-specific Scores")
  
  # function for indexing the multidimensional orthonormal basis 
  eig_index <- function(x, y, z) {
    ind = (x - 1) + HPCA[[4]][[z]][[2]]$FVE * (y - 1) + 1
    return(ind)
  }
  
  # form vectorized versions of the multidimensional orthonormal basis 
  prod_surf <- matrix(list(), nrow = d_tot, 1)
  names(prod_surf) <- groups
  for (g in groups) { 
    # define dimensions of group-specific matrix 
    # prod_surf[[g]] <- data.frame(matrix(
    #   NA, 
    #   nrow = r_tot * f_tot,
    #   ncol = HPCA[[3]][[g]][[2]]$FVE * HPCA[[4]][[g]][[2]]$FVE + 3
    # ))
    
    prod_surf[[g]] <- data.frame(
      reg_index = rep(1:r_tot, each = f_tot), 
      func_index = rep(1:f_tot, r_tot)
    ) %>% 
      mutate(func_reg_index = (reg_index - 1) + 
               r_tot * (func_index - 1) + 1) %>%
      select(-reg_index, -func_index)
    
    for (k in 1:HPCA[[4]][[g]][[2]]$FVE) {
      for (l in 1:HPCA[[3]][[g]][[2]]$FVE) {
        # kth region marginal eigenvector
        v_k <- HPCA[[4]][[g]][[2]]$vectors[, k]
        
        # lth functional marginal eigenfunction
        phi_l <- HPCA[[3]][[g]][[2]]$vectors[, l]
        
        varname <- paste0("X", 3 + eig_index(k, l, g))
        
        # Multdimensional orthonormal basis
        prod_surf[[g]] <- mutate(prod_surf[[g]], 
                                 !!varname := kronecker(v_k, phi_l)) %>% 
          setNames(c("index", map_chr(1:(ncol(.) - 1), ~ paste0("x", .x))))
      }
    }
  }
  
  # create group-specific design matrix 
  data_lmm <- matrix(list(), nrow = d_tot, 1) 
  names(data_lmm) <- groups
  for (g in groups) { 
    data_lmm[[g]] <- filter(data, Group == g)
    data_lmm[[g]] <- inner_join(data_lmm[[g]], prod_surf[[g]], by = "index") 
  }
  
  # calculate variance components, subject-specific scores, and measurement 
  # error variance by linear mixed effects model for each group 
  for (g in groups) { 
    HPCA[[5]][[g]] <- matrix(list(), nrow = 3, 1) 
    names(HPCA[[5]][[g]]) <- c("subj_score", "tau_kl", "sigma_2")
    
    xnam <- paste0("x", 1:(ncol(prod_surf[[g]]) - 1))
    formula <- as.formula(paste("~0+", 
                                paste(xnam, collapse = "+"), 
                                "|Subject")) 
    
    # specify linear mixed effects model 
    RE_struct <- reStruct(formula, pdClass = "pdDiag", data = data_lmm[[g]]) 
    RE_model <- lme(ERP_centered2 ~ - 1, 
                    random = RE_struct, 
                    data = data_lmm[[g]], 
                    control = reControl) 
    
    # extract subject-specific scores 
    HPCA[[5]][[g]][["subj_score"]] <- random.effects(RE_model) 
    
    # extract variance components
    HPCA[[5]][[g]][["tau_kl"]] <- getVarCov(RE_model) 
    
    # extract measurement error variance 
    HPCA[[5]][[g]][["sigma_2"]] <- RE_model$sigma^2 
    
  }
  
  tictoc::toc()
  
  
  # Return Results ---------------------------------------------------------
  
  HPCA[[6]] <- data_lmm
  
  tictoc::toc()
  
  return(HPCA)
}

Code for the 2D HPCA Boostrap procedure:

# install and load required packages
if (!require("pacman")) 
  install.packages("pacman", repos = 'http://cran.us.r-project.org')
p_load(
  "tidyverse", # many useful packages :) 
  "here",      # file management 
  "caTools",
  "refund"
)



HPCA_2d_bootstrap <- function(
  B,       # number of bootstrap samples (scalar)
  func,    # functional domain grid (vector)
  reg,     # regional domain grid (vector)
  HPCA,    # HPCA output from HPCA_decomp.R
  region,  # region to consider for hypothesis test (scalar)
  D,       # list of groups to include in bootstrap (vector)
  N        # group sample sizes (vector)
) {
  
  # Indexing function
  eig_index <- function(x) {
    ind <- (x[1] - 1) + K[d] * (x[2] - 1) + 1
    return(ind)
  }
  
  # 0. Setup Parametric Bootstrap Components -----------------------------
  
  # Define global variables
  f_tot <- length(func)  # Total grid points in functional domain
  r_tot <- length(reg)   # Total grid points in regional domain
  
  # Load fixed effects from HPCA decomposition results
  mu  <- HPCA[[1]]
  eta <- HPCA[[2]]
  
  # Obtain region shift under null hypothesis
  eta_null <- matrix(list(), nrow = length(D))
  
  
  for(d in D) {
    eta_null[[d]] <- matrix(list(), nrow = r_tot)
    for (r in reg) { 
      eta_null[[d]][[r]] <- matrix(0, nrow = f_tot)
    }
  }
  
  
  
  for (d in D) {
    for (r in reg) {
      if (r %in% region) {
        for (d1 in D) {
          eta_null[[d]][[r]] <-
            filter(eta, Electrode == reg[r], Group == d1) %>%
            pull(Eta_Mean) / length(D) + eta_null[[d]][[r]]
        }
      } else {
        eta_null[[d]][[r]] <-
          filter(eta, Electrode == reg[r], Group == d) %>% pull(Eta_Mean)
      }
    }
  }
  
  # Load eigencomponents from HPCA decomposition results
  phi <- matrix(list(), nrow = length(D))
  nu  <- matrix(list(), nrow = length(D))
  
  for(d in D){
    nu[[d]]  <- HPCA[[4]][[d]][[2]]$vectors
    rownames(nu[[d]]) <- reg
    phi[[d]] <- HPCA[[3]][[d]][[2]]$vectors
  }
  
  # Extract number of eigencomponents to include
  K <- rep(NA, length(D)) %>% setNames(D)   
  L <- rep(NA, length(D)) %>% setNames(D)   
  
  for(d in D) {
    K[d] <- HPCA[[4]][[d]][[2]]$FVE
    L[d] <- HPCA[[3]][[d]][[2]]$FVE
  }
  
  # Random effects error variance
  lambda <- matrix(list(), nrow = length(D))
  for(d in D) {
    lambda[[d]] <- diag(HPCA[[5]][[d]][[2]])
  }  
  
  # Measurement error variance
  sigma <- rep(NA, length(D)) %>% setNames(D)
  for(d in D){
    sigma[d] <- HPCA[[5]][[d]][[3]]
  }  
  
  
  # 1. Bootstrap Procedure -----------------------------------------------
  
  # Create list to store group-region shifts for each bootstrap iteration
  eta_b <- matrix(list(), nrow = B)
  for (b in 1:B){
    eta_b[[b]] <- matrix(list(), nrow = length(D))
  }
  
  # Create list to store region shifts under the null for each bootstrap iteration
  eta_b_null = matrix(list(), nrow = B)
  
  # Parametric bootstrap procedure  
  for(b in 1:B) {
    
    # Create storage for group-region specific shifts and null
    for(d in D){
      eta_b[[b]][[d]] = matrix(list(), nrow = r_tot)
    }
    
    eta_b_null[[b]] = matrix(list(), nrow = r_tot)
    
    # Generate data
    y <- matrix(list(), nrow = length(D))
    for (d in D){
      y[[d]] <- matrix(list(), nrow = N[d])
      for (i in 1:N[d]) {
        y[[d]][[i]] <- matrix(list(), nrow = r_tot)
        for (r in reg) {
          noise <- matrix(
            rnorm(f_tot, 0, sqrt(sigma[[d]])), 
            nrow = f_tot
          )
          y[[d]][[i]][[r]] <- mu + eta_null[[d]][[r]] + noise
          for(k in 1:K[d]){
            for(l in 1:L[d]){
              RE_sim <- rnorm(1, 0, sqrt(lambda[[d]][eig_index(c(k, l))]))
              y[[d]][[i]][[r]] <- y[[d]][[i]][[r]] + 
                RE_sim * nu[[d]][r, k] * matrix(phi[[d]][, l])
            }
          }
        }
      }
    }
    
    
    # Estimate mean function
    obs_overall <- matrix(list(), nrow = length(D))
    mu_bs <- map_dfc(
      D, 
      function(g) {
        dat <- map_dfr(
          y[[d]], 
          ~ map_dfc(.x, identity) %>% 
            setNames(reg) %>% 
            mutate(Time = functional) %>% 
            gather(Electrode, ERP, -Time)
        )
        spline_fit <- smooth.spline(dat$Time, dat$ERP)
        group_mean <- predict(spline_fit, functional)$y
      }
    ) %>% 
      # calculate overall mean function as the pointwise average of the group
      # means
      mutate(overall = rowMeans(.)) %>% 
      pull(overall) 
    
    
    eta_b[[b]] <- map(
      D, 
      function(d) { 
        map(
          regional, 
          function(r) {
            dat <- map_dfr(
              y[[d]], 
              ~ map_dfc(.x, ~.x - mu_bs) %>% 
                setNames(reg) %>% 
                mutate(Time = functional) %>% 
                gather(Electrode, ERP_centered, -Time)
            ) %>% 
              filter(Electrode == r)
            spline_fit <- smooth.spline(dat$Time, dat$ERP_centered)
            group_region_mean <- predict(spline_fit, functional)$y
          }
        ) 
      }
    ) %>% 
      setNames(D)
    
    
    # Calculate region shift under the null hypothesis
    for(r in region) {
      eta_b_null[[b]][[r]] <- matrix(0, nrow = f_tot)
      for (d in D) {
        eta_b_null[[b]][[r]] <- eta_b[[b]][[d]][[r]] / length(D) + eta_b_null[[b]][[r]]
      }
    }
    
    #end bootstrap
    print(paste("Bootstrap...", as.character(b)))
    
  }
  
  
  # 3. Calculate p-values ------------------------------------------------
  
  # Find distribution of test statistic under the null hypothesis
  dist <- matrix(list(), nrow = B)
  for (b in 1:B) {
    dist[[b]] <- matrix(0, nrow = r_tot) 
    rownames(dist[[b]]) <- reg
    for (r in region) {
      for (d in D) {
        dist[[b]][r,] <- trapz(func, (eta_b[[b]][[d]][[r]] - eta_b_null[[b]][[r]]) ^ 2) +
          dist[[b]][r,]
      }
    }
    dist[[b]] <- sqrt(dist[[b]])
  }
  
  # Calculate test statistic for observed data
  ts <- matrix(0, nrow = r_tot)
  rownames(ts) <- reg
  for (r in region) {
    for (d in D) {
      eta_data <- HPCA[[2]] %>% filter(Group == d, Electrode == r) %>% pull(Eta_Mean) 
      ts[r, ] <- trapz(func, (eta_data - eta_null[[d]][[r]]) ^ 2) + ts[r, ]
    }
  }
  ts <- sqrt(ts)
  
  # Obtain p-value
  dist_mat <- matrix(unlist(dist), nrow = length(reg))
  rownames(dist_mat) <- reg
  p_val <- vector("logical", length = r_tot)
  names(p_val) <- reg
  for(r in region){
    p_val[r] <- sum(dist_mat[r, ] > ts[r, ]) / B
  }
  
  # 4. Output results ----------------------------------------------------
  
  p_boot <- p_val[region]
  boot <- list(p_boot, eta_b) # return a list with the p-value and common region shift under
  # the null from each run of the bootstrap procedure
  return(boot)
  
}